Over 1,000 Windows 95 Programs

home *** CD-ROM | disk | FTP | other *** search

/ Over 1,000 Windows 95 Programs / Over 1000 Windows 95 Programs (Microforum) (Disc 1).iso / 1249 / root.t < prev next >

Wrap

Text File | 1997-04-18 | 32.9 KB | 1,681 lines

% % "root.t" demonstates techniques for the numerical % solution of non-linear equations % % Sample program for the T Interpreter by: % % Stephen R. Schmitt % 962 Depot Road % Boxborough, MA 01719 % const MATVECT_EPSILON : real := 1.0E-15 const TOLERANCE : real := 0.00000001 const MAX_ITER : int := 50 const DIM : int := 10 type rvector : array[DIM] of real type ivector : array[DIM] of int type rmatrix : array[DIM,DIM] of real type complex : record re : real im : real end record type cvector : array[DIM] of complex program do_bisection do_newton_approx do_newton_exact do_richmond_approx do_richmond_exact do_combined do_brent do_deflation do_lin_bairstow do_newton_two_fncs do_newton_multi_roots do_newton_snle end program %--------------------------------------------------------------------- % demonstrations of the numerical solution of non-linear equations %--------------------------------------------------------------------- % % Fx returns the value of: f(x) = exp(x) - 3*x^2 % % d_Fx returns the first derivative: f'(x) = exp(x) - 6*x % % dd_Fx returns the second derivative: f"(x) = exp(x) - 6 % function Fx( x : real ) : real return exp(x) - 3*x^2 end function function d_Fx( x : real ) : real return exp(x) - 6*x end function function dd_Fx( x : real ) : real return exp(x) - 6 end function const Fx_msg : string := "solving: f(x) = exp(x) - 3*x^2" % % returns value of either of two sumultaneous non-linear equations: % % f0(x0,x1) = x1*exp(x0) - 2 % f1(x0,x1) = x0^2 + x1 - 4 % function Snle( var x : rvector, index : int ) : real case index of value 0: return x[1] * exp( x[0] ) - 2 value 1: return x[0]^2 + x[1] - 4 value: return 1 end case end function const f0_msg : string := "f0(x0,x1) = x1*exp(x0) - 2" const f1_msg : string := "f1(x0,x1) = x0^2 + x1 - 4" % % F1xy returns the value of: f1(x,y) = x^2 + y^2 - 1 % % F2xy returns the value of: f2(x,y) = x^2 - y^2 + 0.5 % function F1xy( x, y : real ) : real return x^2 + y^2 - 1 end function function F2xy( x, y : real ) : real return x^2 - y^2 + 0.5 end function const F1xy_msg : string := "f1(x,y) = x^2 + y^2 - 1" const F2xy_msg : string := "f2(x,y) = x^2 - y^2 + 0.5" % % test Bisection % procedure do_bisection var low, high, root : real low := 2 high := 4 put "bisection method\n" put Fx_msg put "guess interval is: ", low, "...", high if Bisection( low, high, root ) then put "root = ", root else put "could not solve" end if end procedure % % test Newton_approx % procedure do_newton_approx var root : real := 3 var num_roots, i : int put "\nNewton's method for an approximate derivative\n" put Fx_msg put "initial guess is: ", root if Newton_approx( root ) then put "root = ", root else put "could not solve" end if end procedure % % test Newton_exact % procedure do_newton_exact var root : real := 3 put "\nNewton's method for a user provided derivative\n" put Fx_msg put "initial guess is: ", root if Newton_exact( root ) then put "root = ", root else put "could not solve" end if end procedure % % test Richmond_approx % procedure do_richmond_approx var root : real := 3 put "\nRichmond's method with an approximate derivative\n" put Fx_msg put "initial guess is: ", root if Richmond_approx( root ) then put "root = ", root else put "could not solve" end if end procedure % % test Richmond_exact % procedure do_richmond_exact var root : real := 3 put "\nRichmond's method with a user provided derivative\n" put Fx_msg put "initial guess is: ", root if Richmond_exact( root ) then put "root = ", root else put "could not solve" end if end procedure % % test Combined % procedure do_combined var low, high, root : real low := 2 high := 4 put "\nCombined method\n" put Fx_msg put "guess interval is: ", low, "...", high if Combined( low, high, root ) then put "root = ", root else put "could not solve" end if end procedure % % test Brent % procedure do_brent var low, high, root : real low := 2 high := 4 put "\nBrent's method\n" put Fx_msg put "guess interval is: ", low, "...", high if Brent( low, high, root ) then put "root = ", root else put "could not solve" end if end procedure % % test Deflate_polynomial % procedure do_deflation const SMALL : real := 0.00000001 var coeff, roots : rvector var poly_order, num_roots, i : int poly_order := 3 coeff[0] := -40 coeff[1] := 62 coeff[2] := -23 coeff[3] := 1 put "\nDeflating Polynomial method\n" put "the polynomial:" for decreasing i := poly_order...0 do if i > 1 then put coeff[i], "*X^", i, " + "... elsif i = 1 then put coeff[i], "*X + "... elsif i = 0 then put coeff[i]... end if end for put " = 0" if Deflate_polynomial( coeff, 1.1, roots, num_roots, poly_order ) then put "has roots:" for i := 0...poly_order-1 do put "root ", i + 1, " = ", roots[i] end for else put "could not solve" end if end procedure % % test Lin_Bairstow % procedure do_lin_bairstow var coeff : rvector var roots : cvector var r : complex var num_roots, poly_order, i : int poly_order := 3 coeff[0] := 8 coeff[1] := 4 coeff[2] := 2 coeff[3] := 1 put "\nLin-Bairstow method\n" put "the polynomial:" for decreasing i := poly_order...0 do if i > 1 then put coeff[i], "*X^", i, " + "... elsif i = 1 then put coeff[i], "*X + "... elsif i = 0 then put coeff[i]... end if end for put " = 0" if Lin_Bairstow( coeff, roots, poly_order ) then put "has roots:" for i := 0...poly_order-1 do r.re := roots[i].re r.im := roots[i].im put cstr( r ) end for else put "could not solve" end if end procedure % % test Newton2functions % procedure do_newton_two_fncs var rootX, rootY : real rootX := 1 rootY := 3 put "\nNewton's method for two equations\n" put "solving: " put F1xy_msg put F2xy_msg put "initial guess for x is ", rootX put "initial guess for y is ", rootY if Newton2functions( rootX, rootY ) then put "root X = ", rootX put "root Y = ", rootY else put "could not solve" end if end procedure % % test NewtonMultiRoots % procedure do_newton_multi_roots var roots : rvector var num_roots, i : int roots[0] := 2 put "\nNewton's method for multiple roots\n" put Fx_msg put "initial guess is: ", roots[0] if NewtonMultiRoots( roots, num_roots, DIM ) then for i := 0...num_roots-1 do put "root ", i + 1, " = ", roots[i] end for else put "could not solve" end if end procedure % % test NewtonSimNLE % procedure do_newton_snle var Xguess : rvector var num_eqns, i : int num_eqns := 2 Xguess[0] := -0.16 Xguess[1] := 2.7 put "\nNewton's method for simultaneous non-linear equations\n" put "solving for:" put f0_msg put f1_msg for i := 0...num_eqns-1 do put "initial guess for X", i, " = ", Xguess[i] end for if NewtonSimNLE( Xguess, num_eqns ) then put "roots are:" for i := 0...num_eqns-1 do put "X", i, " = ", Xguess[i] end for else put "could not solve" end if end procedure %--------------------------------------------------------------------- % functions for numerical solution of non-linear equations %--------------------------------------------------------------------- % % Solve for a root of a function of a single variable by selecting % an initial interval that contains the root and searching % function Bisection( low, high : real, var root : real ) : boolean if Fx( low ) * Fx( high ) > 0 then return false end if loop exit when abs( high - low ) < TOLERANCE % update guess root := ( low + high ) / 2 if Fx( root ) * Fx( high ) > 0 then high := root else low := root end if end loop return true end function % % Solve for a root of a function of a single variable by selecting % an initial guess and using the function value and an approximate % derivative to search for the solution % function Newton_approx( var root : real ) : boolean var iter : int := 0 var h, diff : real loop h := 0.01 * root if abs( root ) < 1 then h := 0.01 end if % calculate guess refinement diff := 2 * h * Fx( root ) diff := diff / ( Fx( root+h ) - Fx( root-h ) ) % update guess root := root - diff iter := iter + 1 exit when iter > MAX_ITER or abs( diff ) < TOLERANCE end loop if abs( diff ) <= TOLERANCE then return true else return false end if end function % % Solve for a root of a function of a single variable by selecting % an initial guess and using the function value and the % derivative to search for the solution % function Newton_exact( var root : real ) : boolean var iter : int := 0 var diff : real loop % calculate guess refinement diff := Fx( root ) / d_Fx( root ) % update guess root := root - diff iter := iter + 1 exit when iter > MAX_ITER or abs( diff ) < TOLERANCE end loop if abs( diff ) <= TOLERANCE then return true else return false end if end function % % Solve for a root of a function of a single variable by selecting % an initial guess and using the function value and an approximation % of the first and second derivative to search for the solution % function Richmond_approx( var root : real ) : boolean var iter : int var h, diff : real var f1, f2, f3 : real var fd1, fd2 : real loop h := 0.01 * root if abs( root ) < 1 then h := 0.01 end if f1 := Fx( root - h ) % f(x-h) f2 := Fx( root ) % f(x) f3 := Fx( root + h ) % f(x+h) fd1 := ( f3 - f1 ) / ( 2 * h ) % f'(x) fd2 := ( f3 - 2 * f2 + f1 ) / sqrt( h ) % f"(x) % calculate guess refinement diff := f1 * fd1 / ( fd1 ^ 2 - 0.5 * f1 * fd2 ) % update guess root := root - diff iter := iter + 1 exit when iter > MAX_ITER or abs( diff ) < TOLERANCE end loop if abs( diff ) <= TOLERANCE then return true else return false end if end function % % Solve for a root of a function of a single variable by selecting % an initial guess and using the function value and an approximation % of the first and second derivative to search for the solution % function Richmond_exact( var root : real ) : boolean var iter : int var diff, f1 : real var fd1, fd2 : real iter := 0 loop f1 := Fx( root ) fd1 := d_Fx( root ) fd2 := dd_Fx( root ) % calculate guess refinement diff := f1 * fd1 / ( fd1^2 - 0.5 * f1 * fd2 ) % update guess root := root - diff iter := iter + 1 exit when iter > MAX_ITER or abs( diff ) < TOLERANCE end loop if abs( diff ) <= TOLERANCE then return true else return false end if end function % % Brent's method uses interpolation in a variation of % the bi-section method. % function Brent( low, high : real, var root : real ) : boolean const SMALL : real := 0.0000001 % epsilon const VERY_SMALL : real := 0.0000000001 % near zero var iter : int var a, b, c : real var fa, fb, fc : real var d, e, tol : real var small1, small2, small3 : real var p, q, r : real var s, xm : real iter := 1 a := low b := high c := high fa := Fx( low ) fb := Fx( high ) fc := fb % check that the guesses contain the root if fa * fb > 0 then return false end if % start loop to refine the guess for the root loop exit when iter > MAX_ITER iter := iter + 1 if fb * fc > 0 then c := a fc := fa e := b - a d := e end if if abs(fc) < abs(fb) then a := b b := c c := a fa := fb fb := fc fc := fa end if tol := 2 * SMALL * abs(b) + TOLERANCE / 2 xm := ( c - b ) / 2 if abs( xm ) <= tol or abs( fb ) <= VERY_SMALL then root := b return true end if if abs( e ) >= tol and abs( fa ) > abs( fb ) then % perform the inverse quadratic interpolation s := fb / fa if abs(a - c) <= VERY_SMALL then p := 2 * xm * s q := 1 - s else q := fa / fc r := fb / fc q := s * (2 * xm * q * (q - r) - (b - a) * (r - 1)) q := (q - 1) * (r - 1) * (s - 1) end if % determine if improved guess is inside the range if p > 0 then q := -q end if p := abs( p ) small1 := 3 * xm * q * abs(tol * q) small2 := abs(e * q) if small1 < small2 then small3 := small1 else small3 := small2 end if if 2 * p < small3 then % interpolation successfull e := d d := p / q else % use bisection because the % interpolation did not succeed d := xm e := d end if else % use bisection because the range % is slowly decreasing d := xm e := d end if % copy most recent guess to variable a a := b fa := fb % evaluate improved guess for the root if abs(d) > tol then b := b + d else if xm > 0 then b := b + abs( tol ) else b := b - abs( tol ) end if end if fb := Fx( b ) end loop return false end function % % This process combines the bi-section and newton techniques % function Combined( low, high : real, var root : real ) : boolean var iter : int var h, diff : real iter := 0 loop iter := iter + 1 h := 0.01 * root if abs( root ) < 1 then h := 0.01 end if % calculate guess refinement diff := 2 * h * Fx( root ) / ( Fx( root + h ) - Fx( root - h ) ) root := root - diff % check if Newton's method yields a refined guess % outside the range (low, high) if root < low or root > high then % apply Bisection method for this iteration root := ( low + high ) / 2 if Fx( root ) * Fx( high ) > 0 then high := root else low := root end if end if exit when iter > MAX_ITER or abs( diff ) < TOLERANCE end loop if abs( diff ) <= TOLERANCE then return true else return false end if end function % % The deflating polynomial method combines a polynomial reduction % step with Newton's method. As each (real) root is found, the % order of the polynomial to be solved is reduced. % function Deflate_polynomial( var coeff : rvector, initGuess : real, var roots : rvector, numRoots, polyOrder : int ) : boolean var a b, c : rvector var diff, z, X : real var i, iter, n : int iter := 1 n := polyOrder + 1 X := initGuess numRoots := 0 for i := 0...n-1 do a[i] := coeff[i] end for loop exit when iter > MAX_ITER or n <= 1 iter := iter + 1 z := X b[n-1] := a[n-1] c[n-1] := a[n-1] for decreasing i := n-2...1 do b[i] := a[i] + z * b[i+1] c[i] := b[i] + z * c[i+1] end for b[0] := a[0] + z * b[1] diff := b[0] / c[1] X := X - diff if abs(diff) <= TOLERANCE then iter := 1 % reset iteration counter n := n - 1 roots[numRoots] := X numRoots := numRoots + 1 % update deflated roots for i := 0...n-1 do a[i] := b[i + 1] % get the last root if n = 2 then n := n - 1 roots[numRoots] := -a[0] numRoots := numRoots + 1 end if end for end if end loop return true end function % % The Lin-Bairstow method improves on the deflation method to solve % for the complex roots of the following polynomial: % % y = coeff(0) + coeff(1) X + coeff(2) X^2 +...+ coeff(n) X^n % % parameters: % % coeff must be an array with at least order+1 elements. % roots output array of roots % order order of polynomial % function Lin_Bairstow( var coeff : rvector, var roots : cvector, order : int ) : boolean const SMALL : real := 0.00000001 var a, b, c, d : rvector var alfa1, alfa2 : real var beta1, beta2 : real var delta1, delta2, delta3 : real var i, j, k : int var count : int var n : int var n1 : int var n2 : int n := order n1 := n + 1 n2 := n + 2 % is the coefficient of the highest term zero? if abs( coeff[0] ) < SMALL then return false end if for i := 0...n do a[n1-i] := coeff[i] end for % is highest coeff not close to 1? if abs( a[1]-1 ) > SMALL then % adjust coefficients because a(1) != 1 for i := 2...n1 do a[i] := a[i] / a[1] end for a[1] := 1 end if % initialize root counter count := 0 % % start the main Lin-Bairstow iteration loop % initialize the counter and guesses for the % coefficients of quadratic factor: % % p(x) = x^2 + alfa1 * x + beta1 % loop alfa1 := 1 beta1 := 1 loop b[0] := 0 d[0] := 0 b[1] := 1 d[1] := 1 j := 1 k := 0 for i := 2...n1 do b[i] := a[i] - alfa1 * b[j] - beta1 * b[k] d[i] := b[i] - alfa1 * d[j] - beta1 * d[k] j := j + 1 k := k + 1 end for j := n - 1 k := n - 2 delta1 := d[j]^2 - ( d[n] - b[n] ) * d[k] alfa2 := ( b[n] * d[j] - b[n1] * d[k] ) / delta1 beta2 := ( b[n1] * d[j] - ( d[n] - b[n] ) * b[n] ) / delta1 alfa1 := alfa1 + alfa2 beta1 := beta1 + beta2 exit when abs( alfa2 ) < TOLERANCE and abs( beta2 ) < TOLERANCE end loop delta1 := alfa1^2 - 4 * beta1 if delta1 < 0 then % imaginary roots delta2 := sqrt( abs( delta1 ) ) / 2 delta3 := -alfa1 / 2 roots[count].re := delta3 roots[count].im := +delta2 roots[count+1].re := delta3 roots[count+1].im := -delta2 else % roots are real delta2 := sqrt( delta1 ) roots[count].re := ( delta2 - alfa1 ) / 2 roots[count].im := 0 roots[count+1].re := ( delta2 + alfa1 ) / ( -2 ) roots[count+1].im := 0 end if % update root counter count := count + 2 % reduce polynomial order n := n - 2 n1 := n1 - 2 n2 := n2 - 2 % for n >= 2 calculate coefficients of % the new polynomial if n >= 2 then for i := 1...n1 do a[i] := b[i] end for end if exit when n < 2 end loop if n = 1 then % obtain last single real root roots[count].re := -b[2] roots[count].im := 0 end if return true end function % % This function uses Newton's method adapted for the solution of % two simultaneous non-linear equations. % function Newton2functions( var rootX, rootY : real ) : boolean var Jacob : real var fx0, fy0 : real var hX, hY : real var diffX, diffY : real var fxy, fxx : real var fyy, fyx : real var X : real var Y : real var iter : int X := rootX Y := rootY iter := 1 loop hX := 0.01 * rootX if abs( rootX ) < 1 then hX := 0.01 end if hY := 0.01 * rootY if abs(rootY) < 1 then hY := 0.01 end if fx0 := F1xy( X, Y ) fy0 := F2xy( X, Y ) fxx := ( F1xy( X + hX, Y ) - F1xy( X - hX, Y ) ) / 2 / hX fyx := ( F2xy( X + hX, Y ) - F2xy( X - hX, Y ) ) / 2 / hX fxy := ( F1xy( X, Y + hY ) - F1xy( X, Y - hY ) ) / 2 / hY fyy := ( F2xy( X, Y + hY ) - F2xy( X, Y - hY ) ) / 2 / hY Jacob := fxx * fyy - fxy * fyx diffX := ( fx0 * fyy - fy0 * fxy ) / Jacob diffY := ( fy0 * fxx - fx0 * fyx ) / Jacob X := X - diffX Y := Y - diffY iter := iter + 1 exit when iter > MAX_ITER exit when abs( diffX ) < TOLERANCE and abs( diffY ) < TOLERANCE end loop rootX := X rootY := Y if abs( diffX ) <= TOLERANCE and abs( diffY ) <= TOLERANCE then return true else return false end if end function % % This function uses Newton's method adapted for the solution of % multiple simultaneous non-linear equations. % function NewtonMultiRoots ( var roots : rvector, var numRoots : int, maxRoots : int ) : boolean var iter, i : int var h, diff : real var f1, f2, f3 : real var root : real numRoots := 0 root := roots[0] loop iter := 0 loop h := 0.01 * root if abs( root ) < 1 then h := 0.01 end if f1 := Fx( root - h ) f2 := Fx( root ) f3 := Fx( root + h ) if numRoots > 0 then for i := 0...numRoots-1 do f1 := f1 / ( root - h - roots[i] ) f2 := f2 / ( root - roots[i] ) f3 := f3 / ( root + h - roots[i] ) end for end if % calculate guess refinement diff := 2 * h * f2 / ( f3 - f1 ) % update guess root := root - diff iter := iter + 1 exit when iter > MAX_ITER or abs( diff ) <= TOLERANCE end loop if abs(diff) <= TOLERANCE then roots[numRoots] := root numRoots := numRoots + 1 if root < 0 then root := 0.95 * root elsif root > 0 then root := 1.05 * root else root := 0.05 end if end if exit when iter > MAX_ITER or numRoots >= maxRoots end loop if numRoots > 0 then return true else return false end if end function % % This function uses another adaptation of Newton's method for the % solution of multiple simultaneous non-linear equations. It uses % matrix arithmetic functions LUDecomp and LUBackSubst. % function NewtonSimNLE( var X : rvector, numEqns : int ) : boolean var Xdash : rvector var Fvector : rvector var index : ivector var Jmat : rmatrix var i, j : int var moreIter : boolean var iter : int var rowSwapFlag : int var h : real var dummy : boolean iter := 0 loop iter := iter + 1 % copy the values of array X into array Xdash for i := 0...numEqns-1 do Xdash[i] := X[i] end for % calculate the array of function values for i := 0...numEqns-1 do Fvector[i] := Snle( X, i ) end for % calculate the Jmat matrix for i := 0...numEqns-1 do for j := 0...numEqns-1 do % calculate increment in variable number j if abs( X[j] ) > 1 then h := 0.01 * X[j] else h := 0.01 end if Xdash[j] := Xdash[j] + h Jmat[i,j] := ( Snle( Xdash, i ) - Fvector[i] ) / h % restore incremented value Xdash[j] := X[j] end for % j end for % i % solve for the guess refinement vector dummy := LUDecomp( Jmat, index, numEqns, rowSwapFlag ) LUBackSubst( Jmat, index, numEqns, Fvector ) % clear the more-iteration flag moreIter := false % update guess and test for convergence for i := 0...numEqns-1 do X[i] := X[i] - Fvector[i] if abs( Fvector[i] ) > TOLERANCE then moreIter := true end if end for % check iteration limit if moreIter then if iter > MAX_ITER then moreIter := false else moreIter := true end if end if exit when not moreIter end loop return not moreIter end function % % This function performs LU decomposition of a matrix of real values. % function LUDecomp( var A : rmatrix, var Index : ivector, numRows : int, rowSwapFlag : int ) : boolean var i, j : int var k, iMax : int var large, sum : real var z, z2 : real var scaleVect : rvector % initialize row interchange flag rowSwapFlag := 1 % loop to obtain the scaling element for i := 0...numRows-1 do large := 0 for j := 0...numRows-1 do z2 := abs( A[i,j] ) if z2 > large then large := z2 end if end for % j % no non-zero large value? then exit with an error code if large = 0 then return false end if scaleVect[i] := 1 / large end for % i for j := 0...numRows-1 do for i := 0...j-1 do sum := A[i,j] for k := 0...i-1 do sum := sum - A[i,k] * A[k,j] end for A[i,j] := sum end for large := 0 for i := j...numRows-1 do sum := A[i,j] for k := 0...j-1 do sum := sum - A[i,k] * A[k,j] end for A[i,j] := sum z := scaleVect[i] * abs( sum ) if z >= large then large := z iMax := i end if end for % i if j ~= iMax then for k := 0...numRows-1 do z := A[iMax,k] A[iMax,k] := A[j,k] A[j,k] := z end for % k rowSwapFlag := -1 * rowSwapFlag scaleVect[iMax] := scaleVect[j] end if Index[j] := iMax if A[j,j] = 0 then A[j,j] := MATVECT_EPSILON end if if j ~= numRows then z := 1 / A[j,j] for i := j + 1...numRows-1 do A[i,j] := A[i,j] * z end for % i end if end for % j return true end function % % This function uses the LU decomposition of a matrix to solve: % % A x = b % % the solution for x is returned in b. % procedure LUBackSubst( var A : rmatrix, var Index : ivector, numRows : int, var b : rvector ) var i, j : int var idx, k : int var sum : real k := -1 for i := 0...numRows-1 do idx := Index[i] sum := b[idx] b[idx] := b[i] if k > -1 then for j := k...i-1 do sum := sum - A[i,j] * b[j] end for % j elsif sum ~= 0 then k := i end if b[i] := sum end for % i for decreasing i := numRows-1...0 do sum := b[i] for j := i + 1...numRows-1 do sum := sum - A[i,j] * b[j] end for % j b[i] := sum / A[i,i] end for % i end procedure % % returns the absolute value of the argument % function abs( x : real ) : real if x < 0.0 then x := -x end if return x end function % % returns a complex number as a string % function cstr( var c : complex ) : string var s : string var re, im : real re := c.re im := c.im if im >= 0.0 then s := frealstr(re,14,9 ) & " +" & frealstr(im,12,9 ) & "j" else im := -im s := frealstr(re,14,9 ) & " -" & frealstr(im,12,9 ) & "j" end if return s end function